library(Seurat)
library(cowplot)
library(umap)
library(dplyr)
library(Matrix)
library(readxl)
library(openxlsx)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(sctransform)
library(knitr)
library(kableExtra)
#library(biomaRt)
library(DESeq2)
library(escape)
library(dittoSeq)
library(GSEABase)
library(scater)
library(ComplexHeatmap)
gene.lists <- read_xlsx("Munoz_Yilmaz_CellCycle_signatures.xlsx")
gene.lists.names <- colnames(gene.lists)
GOI.lists <- c()
for (i in gene.lists.names) {
tmpList <- gene.lists %>% dplyr::select(all_of(i))
tmpList <- tmpList[!is.na(tmpList)]
GOI.lists[[i]] <- tmpList
}
the initial processing was done with r 3.6.1 with Seurat 3.2.0 - the UMAP comes out slightly differently in r 4.0.3 with Seurat 3.2.3*
#al.dat <- Read10X("200218Yil_data/al/filtered_feature_bc_matrix/")
#f.dat <- Read10X("200218Yil_data/f/filtered_feature_bc_matrix/")
#rf.dat <- Read10X("200218Yil_data/rf/filtered_feature_bc_matrix/")
#rf.rapa.dat <- Read10X("200218Yil_data/rf.rapa/filtered_feature_bc_matrix/")
#al <- CreateSeuratObject(counts = al.dat, project = "al", min.cells = 3, min.features = 200)
#f <- CreateSeuratObject(counts = f.dat, project = "f", min.cells = 3, min.features = 200)
#rf <- CreateSeuratObject(counts = rf.dat, project = "rf", min.cells = 3, min.features = 200)
#rf.rapa <- CreateSeuratObject(counts = rf.rapa.dat, project = "rf.rapa", min.cells = 3, min.features = 200)
#al[["percent.mt"]] <- PercentageFeatureSet(al, pattern = "^mt-")
#f[["percent.mt"]] <- PercentageFeatureSet(f, pattern = "^mt-")
#rf[["percent.mt"]] <- PercentageFeatureSet(rf, pattern = "^mt-")
#rf.rapa[["percent.mt"]] <- PercentageFeatureSet(rf.rapa, pattern = "^mt-")
#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(f, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(rf, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#merged_Seu <- merge(al, c(f,rf,rf.rapa), project = "Diet")
#merged_Seu <- NormalizeData(merged_Seu, normalization.method = "LogNormalize", scale.factor = 10000)
#merged_Seu <- FindVariableFeatures(merged_Seu, selection.method = "vst", nfeatures = 2000)
#merged_Seu <- ScaleData(merged_Seu)
#merged_Seu <- RunPCA(merged_Seu, features = VariableFeatures(object = merged_Seu))
#merged_Seu <- RunUMAP(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- RunTSNE(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- FindNeighbors(merged_Seu, dims = 1:30, verbose = FALSE)
#merged_Seu <- FindClusters(merged_Seu, verbose = FALSE)
#DimPlot(merged_Seu,reduction="umap",group.by="orig.ident",label=TRUE,repel=FALSE)
#FeaturePlot(merged_Seu,reduction="umap",features="mt-Cytb",min.cutoff=0,max.cutoff=4,cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="percent.mt",cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="nFeature_RNA",cols=c("grey","red"))
#save(merged_Seu, file="merged.Robj")
#load("merged.Robj")
#di <- SplitObject(merged_Seu, split.by = "orig.ident")
#for (i in 1:length(di)) {
# di[[i]] <- NormalizeData(di[[i]], verbose = FALSE)
# di[[i]] <- FindVariableFeatures(di[[i]], selection.method = "vst", nfeatures = 2000,
# verbose = FALSE)
#}
#dat.anchors <- FindIntegrationAnchors(object.list=di,dims=1:30)
#integrated <- IntegrateData(anchorset=dat.anchors,dim=1:30)
#DefaultAssay(integrated)<-"integrated"
#integrated <- ScaleData(integrated)
#integrated <- RunPCA(integrated,npcs=30)
#integrated <- RunUMAP(integrated,reduction="pca",dims=1:30)
#integrated <- RunTSNE(integrated,reduction="pca",dims=1:30)
#integrated <- FindNeighbors(integrated, dims = 1:30, verbose = FALSE)
#integrated <- FindClusters(integrated, verbose = FALSE)
NOTE: loading final object to avoid recalculating ssGSEA data
#load("integrated_orig.Robj")
load("../repo_data/integrated_091120.Robj")
#load("../repo_data/integrated_final.Robj")
DefaultAssay(integrated)<-"integrated"
DimPlot(integrated,reduction="umap",split.by="orig.ident",group.by="orig.ident")
DimPlot(integrated,reduction="umap",group.by="orig.ident")
DimPlot(integrated,reduction="umap",group.by="integrated_snn_res.0.8", label=TRUE)
## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
integrated@meta.data$treat_clust <- paste0(integrated@meta.data$orig.ident,integrated@meta.data$integrated_snn_res.0.8)
integrated@meta.data$clust_treat <- paste0(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
integrated@meta.data$celltype <- integrated@meta.data$integrated_snn_res.0.8
Idents(object = integrated) <- integrated$celltype
integrated <- RenameIdents(object = integrated, '16' = '16_Tuft')
integrated <- RenameIdents(object = integrated, '11' = '11_EC')
integrated <- RenameIdents(object = integrated, '13' = '13_EEC')
integrated <- RenameIdents(object = integrated, '2' = '2_Stem')
integrated <- RenameIdents(object = integrated, '5' = '5_Stem')
integrated <- RenameIdents(object = integrated, '10' = '10_Stem')
integrated <- RenameIdents(object = integrated, '14' = '14_Paneth')
integrated <- RenameIdents(object = integrated, '8' = '8_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '9' = '9_Goblet')
integrated <- RenameIdents(object = integrated, '1' = '1_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '4' = '4_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '15' = '15_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '3' = '3_EC_Progenitor')
integrated <- RenameIdents(object = integrated, '6' = '6_EC_Progenitor')
integrated <- RenameIdents(object = integrated, '0' = '0_Early_TA')
integrated <- RenameIdents(object = integrated, '7' = '7_Early_TA')
integrated <- RenameIdents(object = integrated, '12' = '12_Unknown')
integrated[["clust_celltype"]] <- Idents(object = integrated)
Idents(object = integrated) <- integrated$clust_celltype
celltype_colors <- c("2_Stem"="#117733",
"5_Stem"="#999933",
"10_Stem"="#009E73",
"0_Early_TA"="#E69F00",
"7_Early_TA"="#D55E00",
"6_EC_Progenitor"="#0072B2",
"3_EC_Progenitor"="#56B4E9",
"11_EC"="#88CCEE",
"13_EEC"="#6699CC",
"1_Secretory_Progenitor"="#661100",
"4_Secretory_Progenitor"="#882255",
"8_Secretory_Progenitor"="#CC6677",
"15_Secretory_Progenitor"="#AA4499",
"14_Paneth"="#332288",
"16_Tuft"="#000000",
"9_Goblet"="#F0E442",
"12_Unknown"="#999999")
dp.cb <- DimPlot(integrated,reduction="umap", cols=celltype_colors, label=TRUE, repel=TRUE, pt.size=2, label.size=6) + NoLegend()
dp.cb
#pdf('Fig3b.cbSafe.pdf',width=14, height=10)
#dp.cb
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
plotOrder <- c("5","2","10","0","1","3","4","6","7","8","9","11","12","13","14","15","16")
vln_colors <- c("2"="#117733",
"5"="#999933",
"10"="#009E73",
"0"="#E69F00",
"1"="#661100",
"3"="#56B4E9",
"4"="#882255",
"6"="#0072B2",
"7"="#D55E00",
"8"="#CC6677",
"9"="#F0E442",
"11"="#88CCEE",
"12"="#999999",
"13"="#6699CC",
"14"="#332288",
"15"="#AA4499",
"16"="#000000")
Idents(integrated) <- factor(Idents(integrated), levels= plotOrder)
vl.cb <- VlnPlot(integrated, cols=vln_colors, idents = , features = c("Lgr5"), pt.size = 0.5, slot="data")
vl.cb
#pdf('Fig3c_cbSafe.pdf',width=14, height=8)
#vl.cb
#dev.off()
p.integrated <- table(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
round(prop.table(p.integrated,2),3)
##
## al f rf rf.rapa
## 0 0.169 0.191 0.135 0.178
## 1 0.096 0.101 0.135 0.119
## 2 0.082 0.116 0.104 0.102
## 3 0.099 0.091 0.098 0.101
## 4 0.125 0.078 0.082 0.061
## 5 0.076 0.080 0.084 0.072
## 6 0.076 0.074 0.065 0.077
## 7 0.056 0.040 0.096 0.055
## 8 0.078 0.053 0.034 0.037
## 9 0.028 0.046 0.043 0.060
## 10 0.034 0.032 0.044 0.055
## 11 0.026 0.056 0.033 0.029
## 12 0.039 0.028 0.024 0.032
## 13 0.005 0.006 0.009 0.009
## 14 0.005 0.005 0.006 0.007
## 15 0.002 0.002 0.005 0.005
## 16 0.004 0.002 0.004 0.001
DefaultAssay(integrated)<- "integrated"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
dd <- subset(integrated, idents = c("2", "5", "10"), downsample=100)
topvst <- head(VariableFeatures(dd), 500)
mat <- as.matrix(dd@assays$integrated@scale.data) #as.matrix(subset_dd@assays$integrated@scale.data)
mat <- mat[topvst,]
genes <- c(GOI.lists$Biton_S1_ISC.I, GOI.lists$Biton_S1_ISC.II, GOI.lists$Biton_S1_ISC.III)
labels <- c(rep('Biton ISC I', length(GOI.lists$Biton_S1_ISC.I)),
rep('Biton ISC II', length(GOI.lists$Biton_S1_ISC.II)),
rep('Biton ISC III', length(GOI.lists$Biton_S1_ISC.III)))
idxs <- which(genes %in% rownames(mat))
genes <- genes[idxs]
labels <- labels[idxs]
mat <- mat[genes,]
ht <- Heatmap(mat, column_names_side = 'top', row_names_gp = gpar(fontsize = 9), column_names_gp = gpar(fontsize = 12),
column_title = '', name = 'scaled data', column_labels = rep('', ncol(mat)),
column_split = factor(as.character(dd$integrated_snn_res.0.8), levels = c('5', '2', '10')),
row_split = labels, row_order = genes, #column_order = sort(colnames(mat)),
cluster_column_slices = FALSE,
top_annotation = HeatmapAnnotation(cluster = as.character(dd$integrated_snn_res.0.8)))
draw(ht)
#pdf('Fig3D.pdf',width=12, height=10)
#draw(ht)
#dev.off()
data <- as.data.frame(read_excel('cluster_2_5_10stem_gsea.xlsx', sheet = 'Sheet2'))
data
## NAME Comparison NES FDR q-val
## 1 BITON_S1_ISC.I f5_v_al5.noNA.NES 1.3215407 0.2682069200
## 2 BITON_S1_ISC.II f5_v_al5.noNA.NES 2.0391748 0.0000000000
## 3 BITON_S1_ISC.III f5_v_al5.noNA.NES -1.1991149 0.3889467000
## 4 BITON_S1_ISC.I rf5_v_al5.noNA.NES 1.7629898 0.0019750001
## 5 BITON_S1_ISC.II rf5_v_al5.noNA.NES 2.2566988 0.0000000000
## 6 BITON_S1_ISC.III rf5_v_al5.noNA.NES 0.8156711 0.9565237000
## 7 BITON_S1_ISC.I rf.rapa5_v_al5.noNA.NES 0.9629948 0.9024616000
## 8 BITON_S1_ISC.II rf.rapa5_v_al5.noNA.NES 2.0106223 0.0000000000
## 9 BITON_S1_ISC.III rf.rapa5_v_al5.noNA.NES 0.7934456 0.9918505500
## 10 BITON_S1_ISC.I f2_v_al2.noNA.NES 1.3257033 0.2230794400
## 11 BITON_S1_ISC.II f2_v_al2.noNA.NES 2.0813794 0.0000000000
## 12 BITON_S1_ISC.III f2_v_al2.noNA.NES -0.8689138 0.6766603600
## 13 BITON_S1_ISC.I rf2_v_al2.noNA.NES 1.2952416 0.3336595600
## 14 BITON_S1_ISC.II rf2_v_al2.noNA.NES 2.0858900 0.0000000000
## 15 BITON_S1_ISC.III rf2_v_al2.noNA.NES 0.9813662 0.9102961000
## 16 BITON_S1_ISC.I rf.rapa2_v_al2.noNA.NES -1.9235135 0.0047468350
## 17 BITON_S1_ISC.II rf.rapa2_v_al2.noNA.NES 1.8073893 0.0003529412
## 18 BITON_S1_ISC.III rf.rapa2_v_al2.noNA.NES 1.2005840 0.4414399000
## 19 BITON_S1_ISC.I f10_v_al10.noNA.NES 1.0625067 0.6857741500
## 20 BITON_S1_ISC.II f10_v_al10.noNA.NES 2.1364486 0.0000000000
## 21 BITON_S1_ISC.III f10_v_al10.noNA.NES 0.7806379 1.0000000000
## 22 BITON_S1_ISC.I rf10_v_al10.noNA.NES 0.8879838 0.8637682000
## 23 BITON_S1_ISC.II rf10_v_al10.noNA.NES 2.0704706 0.0000000000
## 24 BITON_S1_ISC.III rf10_v_al10.noNA.NES 0.6868493 0.9426303500
## 25 BITON_S1_ISC.I rf.rapa10_v_al10.noNA.NES 0.9626449 0.9918727000
## 26 BITON_S1_ISC.II rf.rapa10_v_al10.noNA.NES 2.0112374 0.0000000000
## 27 BITON_S1_ISC.III rf.rapa10_v_al10.noNA.NES 0.9090225 0.9317172000
data$FDR <- data$`FDR q-val` + 0.001
data$Comparison <- gsub('\\.noNA\\.NES','',data$Comparison)
comparisons <- c('f5_v_al5', 'rf5_v_al5', 'rf.rapa5_v_al5')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl5 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) +
geom_point() +
scale_color_gradient2(midpoint=0, low="blue", mid="white",
high="red", space ="Lab", limits=c(-3,3))+
scale_size_continuous(range=c(1,6))+
ggtitle('Cluster 5 GSEA') + theme_classic() +
theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl5
#pdf('Fig5F.pdf',width=4, height=4)
#cl5
#dev.off()
comparisons <- c('f2_v_al2', 'rf2_v_al2', 'rf.rapa2_v_al2')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl2 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) +
geom_point() +
scale_color_gradient2(midpoint=0, low="blue", mid="white",
high="red", space ="Lab", limits=c(-3,3))+
scale_size_continuous(range=c(1,6))+
ggtitle('Cluster 2 GSEA') + theme_classic() +
theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl2
#pdf('FigS3Da.pdf',width=4, height=4)
#cl2
#dev.off()
comparisons <- c('f10_v_al10', 'rf10_v_al10', 'rf.rapa10_v_al10')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl10 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) +
geom_point() +
scale_color_gradient2(midpoint=0, low="blue", mid="white",
high="red", space ="Lab", limits=c(-3,3))+
scale_size_continuous(range=c(1,6))+
ggtitle('Cluster 10 GSEA') + theme_classic() +
theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl10
#pdf('FigS3Db.pdf',width=4, height=4)
#cl10
#dev.off()
#there is a glitch in this plot, cl2 loses x axis, legend order is different
figS3D <- ggarrange(cl2,cl10, ncol=2, nrow=1)
figS3D
#pdf('FigS3D_incorrect.pdf',width=8, height=4)
#figS3D
#dev.off()
DefaultAssay(integrated)<-"RNA"
colorScheme <- c("#C1BEBF","#fe6100")
fp.Muc2 <- FeaturePlot(integrated,reduction="umap",features="Muc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Muc2
fp.Tff3 <- FeaturePlot(integrated,reduction="umap",features="Tff3",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Tff3
fp.Lyz1 <- FeaturePlot(integrated,reduction="umap",features="Lyz1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lyz1
fp.Defa30 <- FeaturePlot(integrated,reduction="umap",features="Defa30",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Defa30
fp.Chga <- FeaturePlot(integrated,reduction="umap",features="Chga",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Chga
fp.Reg3a <- FeaturePlot(integrated,reduction="umap",features="Reg3a",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Reg3a
fp.Alpi <- FeaturePlot(integrated,reduction="umap",features="Alpi",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Alpi
fp.Atoh1 <- FeaturePlot(integrated,reduction="umap",features="Atoh1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Atoh1
fp.Lgr5 <- FeaturePlot(integrated,reduction="umap",features="Lgr5",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lgr5
fp.Smoc2 <- FeaturePlot(integrated,reduction="umap",features="Smoc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Smoc2
figS4A <- ggarrange(fp.Muc2,fp.Tff3,fp.Lyz1,fp.Defa30,fp.Chga,fp.Reg3a,fp.Alpi,fp.Atoh1,fp.Lgr5,fp.Smoc2, legend = "right", ncol = 5, nrow = 2)
figS4A
#pdf('FigS4A_cbSafe.pdf',width=20, height=10)
#figS4A
#dev.off()
DefaultAssay(integrated)<-"RNA"
gene.names <- rownames(integrated@assays$RNA@data)
Idents(object = integrated) <- integrated$seurat_clusters
stem.subset <- subset(integrated, idents = c("2","5","10"))
levels(stem.subset) <- c("5","2","10")
vln_colors <- c("2"="#117733",
"5"="#999933",
"10"="#009E73")
s.InData <- intersect(gene.names,GOI.lists$mm.s)
g2m.InData <- intersect(gene.names,GOI.lists$mm.g2m)
stem.subset[["percent.mm.s"]] <- PercentageFeatureSet(stem.subset, features = s.InData)
stem.subset[["percent.mm.g2m"]] <- PercentageFeatureSet(stem.subset, features = g2m.InData)
mm.s <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.s",pt.size = 0.3,slot = "data")
mm.s
mm.g2m <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.g2m",pt.size = 0.3,slot = "data")
mm.g2m
figS4b <- ggarrange(mm.s,mm.g2m, legend = FALSE, ncol=2, nrow=1)
figS4b
pdf('figS4b_cbSafe.pdf',width=8, height=4)
figS4b
dev.off()
## png
## 2
integrated[["percent.mm.s"]] <- PercentageFeatureSet(integrated, features = s.InData)
integrated[["percent.mm.g2m"]] <- PercentageFeatureSet(integrated, features = g2m.InData)
mm.g2m.Flag <- if_else(integrated@meta.data$percent.mm.g2m >= 0.3, "Yes", "No")
integrated@meta.data$mm.g2m.Flag <- mm.g2m.Flag
p <- table(integrated$clust_treat,integrated$mm.g2m.Flag)
p.g2m.summary <- round(prop.table(p,2),3)
mm.s.Flag <- if_else(integrated@meta.data$percent.mm.s >= 0.2, "Yes", "No")
integrated@meta.data$mm.s.Flag <- mm.s.Flag
p <- table(integrated$clust_treat,integrated$mm.s.Flag)
p.s.summary <- round(prop.table(p,1),3)
This code is no longer working due to R and package updates but resulting data is stored in the seurat object escape 1.0.0 is probably required but is no longer available. the version in this R image is 1.0.1 escape ssGSEA - run ssGSEA to quantify expression of the BitonI gene set clusters
#egs <- GeneSet(GOI.lists$Biton_S1_ISC.I, setName="BitonI")
#ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
#integrated@meta.data$BitonI_ssGSEA <- ES$BitonI
# egs <- GeneSet(GOI.lists$Biton_S1_ISC.II, setName="BitonII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonII_ssGSEA <- ES$BitonII
# egs <- GeneSet(GOI.lists$Biton_S1_ISC.III, setName="BitonIII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonIII_ssGSEA <- ES$BitonIII
Figure 5G - Biton 1 in cluster5 Figure 5G_II - Biton II in cluster2 Figure 5G_III - Biton III in cluster10 Figure 5H - Pdgfa in cluster5 Figure S3E - Gkn3 in cluster5*
Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust
vln_colors <- c("al5"="#e69f00",
"f5"="#56b4f9",
"rf5"="#117733",
"rf.rapa5"="#d55e00")
plot.order <- c("al5","f5","rf5","rf.rapa5")
Idents(object = clust5.subset) <- factor(Idents(object = clust5.subset), levels = plot.order)
vln.BitonI.ssGSEA <- VlnPlot(clust5.subset, cols = vln_colors, features="BitonI_ssGSEA",slot = "data",pt.size = 0.4) +
stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-I_ssGSEA")
vln.BitonI.ssGSEA
#pdf('Fig5G.pdf',width=4, height=4)
#vln.BitonI.ssGSEA
#dev.off()
vln.Pdgfa <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Pdgfa"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Pdgfa
pdf('FigS4g_cbSafe.pdf',width=4, height=4)
vln.Pdgfa
dev.off()
## png
## 2
vln.Gkn3 <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Gkn3"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Gkn3
pdf('FigS4h_cbSafe.pdf',width=4, height=4)
vln.Gkn3
dev.off()
## png
## 2
clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust
vln_colors <- c("al2"="#e69f00",
"f2"="#56b4f9",
"rf2"="#117733",
"rf.rapa2"="#d55e00")
plot.order <- c("al2","f2","rf2","rf.rapa2")
Idents(object = clust2.subset) <- factor(Idents(object = clust2.subset), levels = plot.order)
vln.BitonII.ssGSEA <- VlnPlot(clust2.subset, cols = vln_colors, features="BitonII_ssGSEA",slot = "data",pt.size = 0.4) +
stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-II_ssGSEA")
vln.BitonII.ssGSEA
#pdf('Fig5G_II.pdf',width=4, height=4)
#vln.BitonII.ssGSEA
#dev.off()
clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust
vln_colors <- c("al10"="#e69f00",
"f10"="#56b4f9",
"rf10"="#117733",
"rf.rapa10"="#d55e00")
plot.order <- c("al10","f10","rf10","rf.rapa10")
Idents(object = clust10.subset) <- factor(Idents(object = clust10.subset), levels = plot.order)
vln.BitonIII.ssGSEA <- VlnPlot(clust10.subset, cols = vln_colors, features="BitonIII_ssGSEA",slot = "data",pt.size = 0.4) +
stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-III_ssGSEA")
vln.BitonIII.ssGSEA
#pdf('Fig5G_III.pdf',width=4, height=4)
#vln.BitonIII.ssGSEA
#dev.off()
bitonPlot <- ggarrange(vln.BitonI.ssGSEA,vln.BitonII.ssGSEA,vln.BitonIII.ssGSEA,nrow=1,ncol=3)
bitonPlot
#pdf('FigS4f_cbSafe.pdf',width=10, height=5)
#bitonPlot
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust
levels(clust5.subset) <- c("al5","f5","rf5","rf.rapa5")
clust5.subset$treat_clust <- factor(x = clust5.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5"))
VlnPlot(clust5.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
vln.Oat.5 <- VlnPlot(clust5.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5
#pdf('Fig5H_Oat.5.pdf',width=4, height=4)
#vln.Oat.5
#dev.off()
clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust
levels(clust2.subset) <- c("al2","f2","rf2","rf.rapa2")
clust2.subset$treat_clust <- factor(x = clust2.subset$treat_clust, levels = c("al2","f2","rf2","rf.rapa2"))
VlnPlot(clust2.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
vln.Oat.2 <- VlnPlot(clust2.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.2
#pdf('Fig5H_Oat.2.pdf',width=4, height=4)
#vln.Oat.2
#dev.off()
clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust
levels(clust10.subset) <- c("al10","f10","rf10","rf.rapa10")
clust10.subset$treat_clust <- factor(x = clust10.subset$treat_clust, levels = c("al10","f10","rf10","rf.rapa10"))
VlnPlot(clust10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
#vln.Oat.10 <- VlnPlot(clust10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
#vln.Oat.10
#pdf('Fig5H_Oat.10.pdf',width=4, height=4)
#vln.Oat.10
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.2.10.subset <- subset(integrated, idents = c("5","2","10"))
Idents(object = clust5.2.10.subset) <- clust5.2.10.subset$treat_clust
levels(clust5.2.10.subset) <- c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10")
clust5.2.10.subset$treat_clust <- factor(x = clust5.2.10.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10"))
VlnPlot(clust5.2.10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
vln.Oat.5.2.10 <- VlnPlot(clust5.2.10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5.2.10
#pdf('Fig5H_Oat.5.2.10.pdf',width=12, height=4)
#vln.Oat.5.2.10
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$orig.ident
levels(integrated) <- c("al","f","rf","rf.rapa")
integrated$orig.ident <- factor(x = integrated$orig.ident, levels = c("al","f","rf","rf.rapa"))
VlnPlot(integrated, features = c("Oat"), split.by = "orig.ident", group.by = "orig.ident", slot= "data",pt.size = 0.5)
vln.Oat.noclust <- VlnPlot(integrated, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.noclust
#pdf('Fig5H_Oat.noclust.pdf',width=4, height=4)
#vln.Oat.noclust
#dev.off()
# DefaultAssay(integrated)<-"RNA"
# Idents(object = integrated) <- integrated$treat_clust
# levels(x = integrated)
# ##
# #adjust the id1 and id2 variables to set up different tests
# #could use a loop for this
# ##
# id1 <- "al5"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
#
# id1 <- "al10"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
#
# id1 <- "al10"
# id2 <- "al5"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
#save(integrated, file="integrated_final.Robj")
sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] ComplexHeatmap_2.6.2 scater_1.18.6
## [3] SingleCellExperiment_1.12.0 GSEABase_1.52.1
## [5] graph_1.68.0 annotate_1.68.0
## [7] XML_3.99-0.5 AnnotationDbi_1.52.0
## [9] dittoSeq_1.2.6 escape_1.0.1
## [11] DESeq2_1.30.1 SummarizedExperiment_1.20.0
## [13] Biobase_2.50.0 MatrixGenerics_1.2.1
## [15] matrixStats_0.58.0 GenomicRanges_1.42.0
## [17] GenomeInfoDb_1.26.7 IRanges_2.24.1
## [19] S4Vectors_0.28.1 BiocGenerics_0.36.1
## [21] kableExtra_1.3.1 knitr_1.31
## [23] sctransform_0.3.2 ggpubr_0.4.0
## [25] ggrepel_0.9.1 ggplot2_3.3.3
## [27] openxlsx_4.2.3 readxl_1.3.1
## [29] Matrix_1.2-18 dplyr_1.0.4
## [31] umap_0.2.7.0 cowplot_1.1.1
## [33] Seurat_3.2.3
##
## loaded via a namespace (and not attached):
## [1] reticulate_1.18 tidyselect_1.1.0
## [3] RSQLite_2.2.3 htmlwidgets_1.5.3
## [5] BiocParallel_1.24.1 Rtsne_0.15
## [7] munsell_0.5.0 codetools_0.2-16
## [9] ica_1.0-2 future_1.21.0
## [11] miniUI_0.1.1.1 withr_3.0.0
## [13] colorspace_2.0-0 highr_0.8
## [15] rstudioapi_0.13 ROCR_1.0-11
## [17] ggsignif_0.6.0 tensor_1.5
## [19] listenv_0.8.0 labeling_0.4.2
## [21] GenomeInfoDbData_1.2.4 polyclip_1.10-0
## [23] farver_2.0.3 pheatmap_1.0.12
## [25] bit64_4.0.5 parallelly_1.23.0
## [27] vctrs_0.6.5 generics_0.1.0
## [29] xfun_0.21 R6_2.5.1
## [31] clue_0.3-58 ggbeeswarm_0.6.0
## [33] rsvd_1.0.3 msigdbr_7.2.1
## [35] locfit_1.5-9.4 bitops_1.0-6
## [37] spatstat.utils_2.0-0 cachem_1.0.3
## [39] DelayedArray_0.16.3 assertthat_0.2.1
## [41] promises_1.1.1 scales_1.1.1
## [43] beeswarm_0.2.3 gtable_0.3.0
## [45] beachmat_2.6.4 Cairo_1.5-12.2
## [47] globals_0.14.0 goftest_1.2-2
## [49] rlang_1.1.4 genefilter_1.72.1
## [51] GlobalOptions_0.1.2 splines_4.0.3
## [53] rstatix_0.6.0 lazyeval_0.2.2
## [55] broom_0.7.4 yaml_2.2.1
## [57] reshape2_1.4.4 abind_1.4-5
## [59] backports_1.2.1 httpuv_1.5.5
## [61] tools_4.0.3 ellipsis_0.3.2
## [63] RColorBrewer_1.1-2 ggridges_0.5.3
## [65] Rcpp_1.0.6 plyr_1.8.6
## [67] sparseMatrixStats_1.2.1 zlibbioc_1.36.0
## [69] purrr_1.0.2 RCurl_1.98-1.2
## [71] rpart_4.1-15 openssl_2.2.0
## [73] deldir_0.2-9 GetoptLong_1.0.5
## [75] viridis_0.5.1 pbapply_1.4-3
## [77] zoo_1.8-8 haven_2.3.1
## [79] cluster_2.1.0 magrittr_2.0.3
## [81] data.table_1.13.6 RSpectra_0.16-0
## [83] scattermore_0.7 circlize_0.4.12
## [85] lmtest_0.9-38 RANN_2.6.1
## [87] fitdistrplus_1.1-3 GSVA_1.38.2
## [89] hms_1.0.0 patchwork_1.1.1
## [91] mime_0.9 evaluate_0.24.0
## [93] xtable_1.8-4 rio_0.5.16
## [95] shape_1.4.5 gridExtra_2.3
## [97] compiler_4.0.3 tibble_3.0.6
## [99] KernSmooth_2.23-17 crayon_1.4.1
## [101] htmltools_0.5.8.1 mgcv_1.8-33
## [103] later_1.1.0.1 tidyr_1.1.2
## [105] geneplotter_1.68.0 DBI_1.1.1
## [107] MASS_7.3-53 car_3.0-10
## [109] cli_3.6.3 igraph_1.2.6
## [111] forcats_0.5.1 pkgconfig_2.0.3
## [113] foreign_0.8-80 scuttle_1.0.4
## [115] plotly_4.9.3 xml2_1.3.2
## [117] vipor_0.4.5 webshot_0.5.2
## [119] XVector_0.30.0 rvest_0.3.6
## [121] stringr_1.4.0 digest_0.6.36
## [123] RcppAnnoy_0.0.18 spatstat.data_2.0-0
## [125] rmarkdown_2.6 cellranger_1.1.0
## [127] leiden_0.3.7 edgeR_3.32.1
## [129] uwot_0.1.10 DelayedMatrixStats_1.12.3
## [131] curl_5.2.1 shiny_1.6.0
## [133] rjson_0.2.20 lifecycle_1.0.4
## [135] nlme_3.1-149 jsonlite_1.8.8
## [137] BiocNeighbors_1.8.2 carData_3.0-4
## [139] limma_3.46.0 viridisLite_0.3.0
## [141] askpass_1.1 pillar_1.4.7
## [143] lattice_0.20-41 fastmap_1.2.0
## [145] httr_1.4.2 survival_3.2-7
## [147] glue_1.4.2 zip_2.1.1
## [149] spatstat_1.64-1 png_0.1-7
## [151] bit_4.0.4 stringi_1.5.3
## [153] blob_1.2.1 BiocSingular_1.6.0
## [155] memoise_2.0.1 irlba_2.3.3
## [157] future.apply_1.7.0
writeLines(capture.output(sessionInfo()), "2024_sessionInfo.txt")